function  [logl,eu,nuhat] = tloglik(para,data) 

global mprior    
  % calculates log-likelihood based on Kalman filter  
  if mprior == 0 
  alpha = para(1);
gamma = para(2);
beta1  = para(3);
kappa1 = para(4);
tau1   = para(5);
rho_a  = para(6);
rho_f  = para(7);
rho_i  = para(8);
sd_a   = para(9);
sd_f   = para(10);
sd_i   = para(11);
gy     = para(12);
lnA_0  = para(13);
yst    = para(14);
lnP_0  = para(15); 
sd_p   = para(16); 

bound  = zeros(length(para),2);
bound(14,1) = -100; 
bound(:,2) = [100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100]; 
crit1  = para-bound(:,1); crit2=bound(:,2)-para; 
                 if min(crit1)<0 || min(crit2)<0
logl = -1000000; 
                 return 
                 else
iss    = gy-log(beta1) ;


A = [1 1/tau1;
     0 beta1]; 

B = [1+gamma/tau1 alpha/gamma;
     -kappa1  1]; 
     
eu = 0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  determinacy condition violated  *************')
    eu = 1;
end    


if eu==0
A_Aa_mat = [1-rho_a, -rho_a/tau1, 1/tau1; 
            -kappa1,  1-beta1*rho_a,  0;
            -gamma, -alpha, 1];      
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;0;0];

A_Af_mat = [1-rho_f, -rho_f/tau1,  1/tau1 ;
            -kappa1,  1-beta1*rho_f, 0;
            -gamma, -alpha,  1];         
solvec_f = inv(A_Af_mat) * [0;kappa1;0];

A_Ai_mat = [1-rho_i,  -(rho_i)/tau1,  1/tau1;
             -kappa1,  1-beta1*rho_i, 0;
             -gamma, -alpha,  1];        
solvec_i = inv(A_Ai_mat) * [0;0;1];
 ya = solvec_a(1);
 pia = solvec_a(2);
 ia = solvec_a(3);
 yf = solvec_f(1);
 pif = solvec_f(2);
 i_f = solvec_f(3);
 yi = solvec_i(1);
 pii = solvec_i(2);
 ii=solvec_i(3);


%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho       = diag([rho_a;rho_f;rho_i]);

sd       = [sd_a;sd_f;sd_i];


% Augment the state vector 
rho = diag([diag(rho);1;1]);
rho(4,1) = rho_a; 
TC  = [0;0;0;gy;0]; 
sd  = [sd;sd_a;sd_p];
Vs = diag([sd_a^2/(1-rho_a^2);sd_f^2/(1-rho_f^2);sd_i^2/(1-rho_i^2);0;0]);


% system matrices for the macro part 
mac = [ya yf yi 1 0; pia pif pii 0 1;ia i_f ii 0 1]; 


% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 

Az      = [ yst ; 0 ; iss]; 
Bz      = mac;  

% condition on the first observation 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 
s00(5)   = lnP_0;
llik_c   = zeros(nobs,1);       % log likelihood contribution
target   = zeros(nobs,1); 
s0t       = zeros(1,ns);    
ltrvecP0 = zeros(1,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 
nu        = data(1,:)'-(Az+Bz*(TC+rho*s00)); 
Ft        = Bz*Vs*Bz'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz'*inv(Ft)*(Vs*Bz')'))';

s1t = s00'; 
nx    = size(s1t,1);   
target(1) = s00(end); 

t=2;
% log likelihood by Kalman filter
while t<nobs+1
     
   % forecasting
   
   alphahat = TC+rho*s1t';
   P1t      = ltrmat((ltrvecP1(1,:))');
   Phat1    = rho*P1t*rho'+(diag(sd))^2; 
   yhat    = Az + Bz*alphahat;
     
   nu      = data(t,:)-yhat'; 
   Ft       = Bz*Phat1*Bz';
 
        
  % updating step 
  s1t   = (alphahat+Phat1*Bz'*inv(Ft)*nu')';
  ltrvecP1 =(ltrvec(Phat1-Phat1*Bz'*inv(Ft)*Bz*Phat1'))'; 
  nuhat(t,:) = data(t,:)-yhat';
  target(t) = s1t(:,end);  
  
  % llik_c computation 
   
      llik_c(t) = -0.5*nz*log(2*pi) - 0.5*log(det(Ft)) - 0.5*nu*inv(Ft)*nu';
   
  t=t+1;
end

logl = sum(llik_c);
elseif eu>0  
logl = -1000000;
 end                  % determinacy case 
end                   %bound 
  end 

  elseif mprior == 1   
alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a  = para(8);
rho_f  = para(9);
rho_i  = para(10);
sd_a   = para(11);
sd_f   = para(12);
sd_i   = para(13);
piess  = para(14);
gy     = para(15);
sigu   = para(16);
lnA_0  = para(17);
yst    = para(18);
p11    = para(19);
p22    = para(20);

bound  = zeros(length(para),2);
bound(17:18,1) = -100; 
bound(:,2) = [100;100;100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;100]; 

crit1  = para-bound(:,1); crit2=bound(:,2)-para; 
                 if min(crit1)<0 || min(crit2)<0
logl = 1000000; 
                 return 
                 else
iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------

A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

A1_fwz = [1 (1/tau1);
           0 1];

A2_fwz = [1 (1/tau1);
          0 1];
      
B1_fwz = [1+(gamma1/tau1) (alpha1/tau1)
         -kappa1*(1/beta1) (1/beta1)]; 

B2_fwz = [1+(gamma2/tau1) (alpha2/tau1)
         -kappa1*(1/beta1) (1/beta1)];      
     
eu = 0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu = 1;
end    


if eu==0
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho       = diag([rho_a;rho_f;rho_i]);
lambda_a1 = tau1*y1a+1+pi1a;
lambda_a2 = tau1*y2a+1+pi2a;
lambda_f1 = tau1*y1f+pi1f;
lambda_f2 = tau1*y2f+pi2f;
lambda_i1 = tau1*y1i+pi1i;
lambda_i2 = tau1*y2i+pi2i;
lambda1   = [lambda_a1;lambda_f1;lambda_i1];
lambda2   = [lambda_a2;lambda_f2;lambda_i2];


if mprior <3  % homoskedasticity    
sd       = [sd_a;sd_f;sd_i];
sd1 = sd; sd2 = sd; 
sd_i1=sd_i; sd_i2=sd_i; 

elseif mprior >=3 % regime-dependent heteroskedasticity 
sd1 = [sd_a;sd_f;sd_i1]; 
sd2 = [sd_a;sd_f;sd_i2];  
end 

A1 = zeros(1,1);
A2 = zeros(1,1);
B1 = zeros(1,3);
B2 = zeros(1,3);

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho = diag([diag(rho);1]);
rho(4,1) = rho_a; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a];
sd2  = [sd2;sd_a];  
Vs1 = diag([sd_a^2/(1-rho_a^2);sd_f^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs2 = diag([sd_a^2/(1-rho_a^2);sd_f^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = [ p11, 1-p11;1-p22, p22];  
trans = trans';
probs = [(1-p22)/(2-p11-p22);(1-p11)/(2-p11-p22)];

Vs  = probs(1)*Vs1+probs(2)*Vs2; 

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 3;              
nxmax = 2^nhist; 

Az      = [ yst ; piess ; iss]; 
Az1 = Az; Az2 = Az; 
Bz      = [ mac1 ; B1(1,:) 0]; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 
Cz     = diag([0;0;sigu*ones(nz-2,1)]); 

% condition on the first observation 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz*(TC+rho*s00)); 
Ft        = Bz*Vs*Bz'+Cz*Cz'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz'*inv(Ft)*(Vs*Bz')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(2*nx,1);      
s0t       = zeros(2*nx,ns);    
ltrvecP0 = zeros(2*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat = TC+rho*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho*P1t*rho'+(diag(sd1))^2; 
   Phat2    = rho*P1t*rho'+(diag(sd2))^2; 
   yhat1    = Az1 + Bz1*alphahat;
   yhat2    = Az2 + Bz2*alphahat; 
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   
   Ft1      = Bz1*Phat1*Bz1'+Cz*Cz';
   Ft2      = Bz2*Phat2*Bz2'+Cz*Cz'; 
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*inv(Ft2)*nu2'; 
   
   % updating step 
   s0t(i,:)    = (alphahat+Phat1*Bz1'*inv(Ft1)*nu1')';
   s0t(nx+i,:) = (alphahat+Phat2*Bz2'*inv(Ft2)*nu2')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*inv(Ft2)*Bz2*Phat2'))'; 
   
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx))]; 
   
   % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (2*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat = kron(eye(nx),[1 1]); 
   sumprobx2 = selmat*probx2;
   sumzeroprob = selmat*zeroprob;
   probscore  = sumprobx2+sumzeroprob;
   [B,index] =sortrows(probscore);
   collind   = index(1:2*nx-nxmax-sum(zeroprob)); 
   
   
   s1t = zeros(nxmax,ns);
   ltrvecP1 = zeros(nxmax,ntr); 
   probx = zeros(nxmax,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*2+1; 
       
       if sumzeroprob(i)>0; 
        % no merge, just eliminate 
        
         if zeroprob(i2)==0 
            s1t(j,:)=s0t(i2,:);
            ltrvecP1(j,:)=ltrvecP0(i2,:);
            probx(j,:) = probx2(i2); 
            j=j+1;
         end 
         
         if zeroprob(i2+1)==0
            s1t(j,:)=s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvecP0(i2+1,:); 
            probx(j,:)=probx2(i2+1);
            j=j+1;
         end 
         
       else
           
         x=i==collind;
         x=sum(x);
         if x==1
         % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) + (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))-s1t(j,:)'*s1t(j,:))';
            j=j+1;
         else
          % don't collapse density 
             s1t(j,:)=s0t(i2,:);
             ltrvecP1(j,:) = ltrvecP0(i2,:); 
             probx(j,:) = probx2(i2); 
             j=j+1; 
             s1t(j,:)=s0t(i2+1,:);
             ltrvecP1(j,:) = ltrvecP0(i2+1,:);
             probx(j,:) = probx2(i2+1);
              j=j+1;
         end

       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
t=t+1;
end

logl = sum(llik_c);
elseif eu>0  
logl = 1000000;
 end                  % determinacy case 
end                   %bound 
   
   elseif mprior==2 
   
   alpha1 = para(1);
gamma1 = para(2);
beta1  = para(3);
kappa1 = para(4);
tau1   = para(5);
rho_a1   = para(6);
rho_f  = para(7);
rho_i  = para(8);
sd_a1   = para(9);
sd_f1   = para(10);
sd_i1   = para(11);
sd_a2   = para(12);
sd_f2   = para(13);
sd_i2   = para(14);
piess  = para(15);
gy     = para(16);
lnA_0  = para(17);
yst    = para(18);
p11    = para(19);
p22    = para(20);

bound  = zeros(length(para),2);
bound(6:8,1)=0.0001;
bound(17:18,1) = -100; 
bound(:,2) = [100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;100;0.9999;0.9999]; 


crit1  = para-bound(:,1); crit2=bound(:,2)-para; 
                 if min(crit1)<0 || min(crit2)<0
logl = -1000000; 
                 return 
                 else
iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------

alpha2 = alpha1; gamma2=gamma1; 

A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

     
eu = 0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu = 1;
end

rho_a2 = rho_a1; 

                   if eu==0
A_Aa_mat = [1-p11*rho_a1, -(1-p11)*rho_a2, -(p11*rho_a1)/tau1, -((1-p11)*rho_a2)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a1, 1-p22*rho_a2, -((1-p22)*rho_a1)/tau1, -(p22*rho_a2)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a1, -beta1*(1-p11)*rho_a2, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a1, 1-beta1*p22*rho_a2, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [(p11*rho_a1+(1-p11)*rho_a2)/tau1;((1-p22)*rho_a1+p22*rho_a2)/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute term structure
%----------------------------------------------------------------------
rho1      = diag([rho_a1;rho_f;rho_i]);
rho2      = diag([rho_a2;rho_f;rho_i]);
lambda_a1 = tau1*y1a+1+pi1a;
lambda_a2 = tau1*y2a+1+pi2a;
lambda_f1 = tau1*y1f+pi1f;
lambda_f2 = tau1*y2f+pi2f;
lambda_i1 = tau1*y1i+pi1i;
lambda_i2 = tau1*y2i+pi2i;
lambda1   = [lambda_a1;lambda_f1;lambda_i1];
lambda2   = [lambda_a2;lambda_f2;lambda_i2];

sd1 = [sd_a1;sd_f1;sd_i1]; 
sd2 = [sd_a2;sd_f2;sd_i2];  
 

A1 = zeros(1,1);
A2 = zeros(1,1);
B1 = zeros(1,3);
B2 = zeros(1,3);

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho1 = diag([diag(rho1);1]);
rho2 = diag([diag(rho2);1]); 
rho1(4,1) = rho_a1;
rho2(4,1) = rho_a2; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a1]; sd2=[sd2;sd_a2]; 

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = [ p11, 1-p11;1-p22, p22];  
trans = trans'; 
probs = [(1-p22)/(2-p11-p22);(1-p11)/(2-p11-p22)];


Vs1 = probs(1)*diag([sd_a1^2/(1-rho_a1^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0])+probs(2)*diag([sd_a2^2/(1-rho_a2^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 2^nhist; 

Az      = [ yst ; piess ; iss];
Az1     = Az;
Az2     = Az;
Bz      = [ mac1 ; B1(1,:) 0]; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 


% condition on the first observation 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz*(TC+rho1*s00)); 
Ft        = Bz*Vs1*Bz'; 
iFt       = Ft\eye(3);
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*iFt*nu;

s00      = s00+(Vs1*Bz')*iFt*nu; 
ltrvecP1 = (ltrvec(Vs1-Vs1*Bz'*iFt*(Vs1*Bz')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(2*nx,1);      
s0t       = zeros(2*nx,ns);    
ltrvecP0 = zeros(2*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat1 = TC+rho1*s1t(i,:)';
   alphahat2 = TC+rho2*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2;
   Phat2    = rho2*P1t*rho2'+(diag(sd2))^2;
   yhat1    = Az1 + Bz1*alphahat1;
   yhat2    = Az2 + Bz2*alphahat2; 
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   iFt1     = Ft1\eye(3);
   iFt2     = Ft2\eye(3); 
   
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*iFt1*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*iFt2*nu2'; 
   
   % updating step 
   s0t(i,:)    = (alphahat1+Phat1*Bz1'*iFt1*nu1')';
   s0t(nx+i,:) = (alphahat1+Phat2*Bz2'*iFt2*nu2')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*iFt1*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*iFt2*Bz2*Phat2'))'; 
   
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx))]; 
   
   % collapse the state vector 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (2*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat = kron(eye(nx),[1 1]); 
   sumprobx2 = selmat*probx2;
   sumzeroprob = selmat*zeroprob;
   probscore  = sumprobx2+sumzeroprob;
   [B,index] =sortrows(probscore);
   collind   = index(1:2*nx-nxmax-sum(zeroprob)); 
   
   
   s1t = zeros(nxmax,ns);
   ltrvecP1 = zeros(nxmax,ntr); 
   probx = zeros(nxmax,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*2+1; 
       
       if sumzeroprob(i)>0; 
        % no merge, just eliminate 
        
         if zeroprob(i2)==0 
            s1t(j,:)=s0t(i2,:);
            ltrvecP1(j,:)=ltrvecP0(i2,:);
            probx(j,:) = probx2(i2); 
            j=j+1;
         end 
         
         if zeroprob(i2+1)==0
            s1t(j,:)=s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvecP0(i2+1,:); 
            probx(j,:)=probx2(i2+1);
            j=j+1;
         end 
         
       else
           
         x=i==collind;
         x=sum(x);
         if x==1
         % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) + (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))-s1t(j,:)'*s1t(j,:))';
            j=j+1;
         else
          % don't collapse density 
             s1t(j,:)=s0t(i2,:);
             ltrvecP1(j,:) = ltrvecP0(i2,:); 
             probx(j,:) = probx2(i2); 
             j=j+1; 
             s1t(j,:)=s0t(i2+1,:);
             ltrvecP1(j,:) = ltrvecP0(i2+1,:);
             probx(j,:) = probx2(i2+1);
              j=j+1;
         end

       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
t=t+1;
end


logl = sum(llik_c);
%logl = llik_c; 
                   elseif eu>0  
logl = -1000000;
 end                  % determinacy case 
end                   %bound 
if imag(logl)==0;
logl = logl;
else 
logl = -1000000;
end  
   
   
   end 